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Abstract 

We present the method of moments approach to pricing barrier-type op- 
tions when the underlying is modelled by a general class of jump diffusions. 
By general principles the option prices are linked to certain infinite dimen- 
sional linear programming problems. Subsequently approximating those sys- 
tems by finite dimensional linear programming problems, upper and lower 
bounds for the prices of such options are found. As numerical illustration 
we apply the method to the valuation of several barrier-type options (dou- 
ble barrier knockout option, American corridor and double no touch) under a 
number of different models, including a case with deterministic interest rates, 
and compare with Monte Carlo simulation results. In all cases we find tight 
bounds with short execution times. Theoretical convergence results are also 
provided. 
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1 Introduction 



Barrier and barrier-type options are among the most widely and frequently traded 
exotic options, especially in the area of Foreign Exchange, which makes their 
valuation an important topic. For example, a double barrier option is cancelled 
depending on whether or not two levels have been crossed before maturity. Since 
the pay-off of a barrier option depends on the entire path of the underlying, it is 
clear that its valuation is more involved than that of a standard European type 
option. 

A well-documented empirical observation is that financial returns data typically 
possess features such as asymmetry, heavy tails and excess kurtosis, which cannot 
be captured by the classical geometric Brownian motion model (GBM). Related 
is the well known fact that under the GBM model it is not possible to calibrate 
option prices to the volatility surface. One of the successful modifications that has 
been proposed is to introduce jumps in the evolution and work with Levy models. 
Popular examples of such Levy models are VG, CGMY, NIG, GH and KoBoL. 
This approach is classical by now and we refer to the standard references [16 , [3j 
and d] for further financial motivations for the use of jump models, background 
and references. In a separate development (see e.g. [l] and [9j) it was noted that 
commodity prices often display features such as mean-reversion and jumps that 
are clearly not captured by the geometric Brownian motion model, and it was 
proposed to employ jump-diffusion models to incorporate those effects. 

The valuation of barrier options has attracted a good deal of attention and there 
exists currently a body of literature dealing with different aspects of pricing barrier 
options. In particular, for double barrier options, and [15] developed a Laplace 
transform approach in the geometric Brownian motion setting. [T7] derived semi- 
analytical expressions in a jump-diffusion setting with exponential jumps, also 
using a transform approach. [2] considered double no touches in a setting with 
exponential jumps, allowing the process dynamics to change after a barrier is 
breached. [7| used eigenfunction expansions to price double barrier options in a 
CEV setting. 

The mentioned papers exploit specific features of the model under consideration 
and can therefore not be readily generalized and applied to a different settings. A 
general approach, based on a characterization of the moments of the underlying 
process, was followed by [13] to price a class of exotic options. In a diffusion setting 
[13j derived upper and lower bounds for the price of exotic options in terms of 
semidefinite programs, and provided theoretical and numerical convergence results 
for these bounds. Before that, using linear programming, [TT] developed a method 
of moments algorithm to calculate first exit time probabilities and moments of a 
diffusion. 

In this paper we will follow a method of moments approach to price double bar- 
rier options in a general setting of a polynomial-type jump-diffusion. We will 
also allow the rate of discounting, which is typically taken to be constant, to be 
a function of time and underlying. We will now briefly describe the method of 
moments approach. The first step is to express the option price as an integral 
with respect to two measures, the (discounted) expected exit measure and the 
(discounted) expected occupation measure. The former describes the law of the 
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underlying at expiration or at crossing the barrier, while the latter described the 
law of the process until this moment. Restricting ourselves to pay-ofFs that are 
piecewise polynomial functions of the underlying, the value can then be expressed 
as a linear combination of moments of these two measures. The moments of those 
two measures are subsequently shown to satisfy an infinite dimensional linear sys- 
tem. To the price can thus be associated the two linear programming problems 
of minimization and maximization of the latter criterion over the spaces of mea- 
sures. By adding conditions on the moments that guarantee that a given sequence 
is equal to the moments of a measure, one is led to an infinite dimensional linear 
programming problem or a semi-definite programming problem. By restricting to 
a finite number of moments we arrive at a finite dimensional linear programming 
problem or a semi-definite programming problem. 

We will numerically illustrate this method for a American corridor and double 
no touch and double knock-out option under different models, by solving linear 
programming problems. In all cases we find tight bounds, with short execution 
times. We also provide a convergence proof to show that the values of the linear 
programming problems converge monotonically to the value of the option if the 
number of moments employed is increased. 

The remainder of the paper is organized as follows. In Section [2] we specify the 
model and the problem setting. Section [3] is devoted to the method of moments, 
describes the algorithm and provides a convergence proof. Section H] provides the 
implementation and numerical examples. Proofs are deferred to the Appendix. 

2 Problem setting 

Assume that the underlying X evolves according to the SDE 



where 1^ is a Brownian motion and J is a pure jump Levy process (that is, a 
process with independent stationary increments without Gaussian component), 
and cr,6, and A are given functions that will be specified below. In this setting 
we will value a barrier option of knock-out type with pay-off h{T,XT) at the 
maturity time T if the underlying has not left a set B before time T and that pays 
a stream of payments g{s,Xs) until the first moment tb that X leaves B or time 
r, whichever comes earlier. Modelling the risk neutral discounting as a function 
r = r{t,Xt) of t and Xf it follows by standard arbitrage pricing principles that 
the value v of this contract is given by 



dXt = b{t, Xt)dt + a{t, Xt)dWt + A(t, Xt)dJu Xo = xq 






where tb = inf{t >0:Xt^B} and 




(3) 
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We will restrict ourselves to the case that the functions h, g are piecewise polyno- 
mial functions, that is, for some partitions {Cj} and {Di} of [0,T] x M, 



h{t,x) = '^hi{t,x)l{(^t,x)ec,}, 9{t,x) = '^9i{t,x)I{(t,x)&D,}, (4) 

i=l 1=1 

where hi,gi are polynomials in {t,x). Note that many contracts have a pay-off 
function that is of this form, including call and put options and straddles. We 
observe that we will then be able to express the value v in terms of moments of 
certain probability measures, reducing the calculation of v to the calculation of 
these moments. Further, we will assume that X is a 'polynomial' process, that is, 
in eqs. ([I]) and ^ 

b{t,x), a'^{t,x), X{t,x) and r{t,x) are polynomials, 

such that ([1]) admits a unique (weak) solution. Associated to X is the infinitesimal 
generator that acts on functions / in its domain as 

where Bf is an integro-differential operator given by 



Bf{t,x) 



df 

fix + X{t, x)y) - fit, x) - X{t, ^)-Q^{t, x)yl{\y\<i} 



A(dy), (6) 



where A denotes the Levy jump measure of X. Note that the operator A maps 
polynomials to polynomials, which is an essential property needed in the moment 
approach, as shown in the next section. 

We next present some models that are included in our setting. 



Examples. 

• The classical geometric Brownian motion satisfies the SDE 

dX = bXdt + aXdW 

where W denotes a one-dimensional Brownian motion, and has the infinites- 
imal generator 

A/(x) = 6x- + ^-^. (7) 

• Levy models (for an overview see e.g. [T6| or [3]) For example, the variance 
gamma process ([H]) evolves according to dX = bidt + dZ, where bi is a 
constant and Z is a Levy process with Levy measure 

rjidx) = ^e-^^-/{,>o}d^ + |^e-^l-ll{,<o}d^, (8) 

with C, G, M positive constants, and its infinitesimal generator is specified 
by 

Afix) = b^ + J [fix + y)- fix)] vidx). (9) 
where b = bi — J^^ xrj{dx). 
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Additive processes with polynomial time-dependent coefficients (see e.g. ^ 
for background), obtained by taking b, cj^ and A to be polynomials of t only. 
For example, X evolving according to dX = bi{t)dt + dZ for a Levy process 
Z and a polynomial bi[t). 

Affine processes, obtained by taking A constant and 6, o"'^ affine functions in 
X, independent of t (see e.g. [5] for applications of affine models in finance). 
An example of an affine diffusion is the Cox Ingersoll Ross (CIR) model, 
which is a mean-reverting diffusion satisfying the SDE 

dX = a{b - X)dt + (jy/XdW, a,b>0, (10) 

with the infinitesimal generator 

Am = a{b-x)- + —^. (11) 



3 Method of moments 



Denoting by v and /i the discounted exit location measure and the discounted 
occupation measure given by 

K^)=i^[e-"^/{(.,x.)eA}], KA)=E 

for Borel sets A £ B{[0,T] x M), the value v of the contract can be expressed as 

= y h{t,x)v{dt,dx) + j g{t,x)fi{dt,dx) (12) 

The measure /i describes the distribution of the process {t, Xt) before the stopping 
time T whereas the measure v describes the distribution upon termination at r. 
For example, in the case of a up-and-out barrier option at level B^, termination 
occurs if the barrier is crossed or the maturity T is reached. 

In view of the form ^ of g and h, v can be expressed in terms of the moments of 
H and u, as follows: 

- = EE^M-M-i? +EE^m(-)/^&\ (13) 

where we denote by mij = J f x^m{dt, dx) the ijth moment of a measure m and 
by = u[- n Cm) and /x^™) = /u(- Pi Dm) the restrictions of u and /i to Cm and 
and where dij{m) and bij{m) are some constants. 



' I{{s,Xs)£A}ds 



3.1 The adjoint equation 



The measures fi and are closely related to each other and to the generator of the 
underlying process X. Informally, for suitably regular / and all bounded stopping 
times T Dynkin's lemma yields that 

E[e-^^fiT,Xr)]-E[f{0,Xo)]=E 



e-^\Af -rf){t,Xt)dt 



(14) 
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where Af is given in ([5]), which can be expressed in terms of the measures v and 
/X as 

^ fit, x)uidt, dx) = /(O, xo) + l{Af- rf){t, x)^i{dt, dx). (15) 



The identity ([T5]) is called the basic adjoint equation (See e.g. [II]). As noted 
before, a formal application of the generator shows that A maps polynomials to 
polynomials. More specifically, by applying (jlSp to a monomial fij{t,x) = fx^ we 
obtain the following infinite system of equations linking the moments of yu and i': 



J fx^uldtjdx) — xlli=Q = ^^^Ck,i{i,j) J t^x'' fi{dt,dx). 



(16) 



k,l 

l^k 



where {Afij — rfij){t,x) = i'^k,i{i',j)t x or equivalently, in compact notation. 



\li=o = y^^Ck,iii,j)fJ-k,i- 



(17) 



kJ 



The following result provides sufficient conditions to justify this informal analysis: 



Proposition 1 Suppose that for all k = 0,1,2, . 



\y\^ilAy^)Aidy) + E 



-at 



X.fdt 



< oo. 



(18) 



Then eqn. (I17p holds for all i, j = 0,1,2, .. .. 
The proof is deferred to the Appendix. 

Remark. Partial barrier or forward starting barrier options can also be included 
in this setting by slightly adapting the definitions. With f = inf{i G [Tq, T] : Xt ^ 
B} it holds that 



E[e~''^fif,Xf)]- E[f{To,XT,)] = E 
which leads to the adjoint equation 



To 



e-''^{Af-af)it,Xt)dt 



I fdi> = I fduo + /(^/ - "Z)^-"- 



3.1.1 Truncation 



We restrict ourselves now to contracts that are knocked out if X leaves a finite 
interval, so that v is given by ^ with h{t,x) = for x ^ B := If the 

minimum A of \{t,x) over [0, T] x is strictly positive, it is possible to 

derive for such double knock out contracts a modification of the adjoint equations 
that is valid without integrability restrictions. To that end, note that a double 
knock-out option becomes worthless at the first time that a jump occurs of size 
larger than L_|_ := (6+ — b-)/\ or smaller than L_ := (6_ — 6+)/A, since any 
such jump will take X out of the interval 6+]. As such jumps occur at a rate 
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A* = A(M\[L_, independent of the smaller size jumps and the diffusion part, 
it follows that the value v of the contract does not change if we replace A and r 

by 

A = A(- n [L_,L+]), f = r + A*, 

which corresponds to replacing the underlying X by the process X that is 'killed' 
when the first jump occurs with size larger than L-|_ or smaller than L_. In 
summary, using , to denote the ijth moments of the exit and occupation 
measures v, fi of the killed process X, we have the following result (with a proof 
in the Appendix): 

Corollary 1 If B = and A > 0, then v is given by (|12p where v^^ and fi"^^ 

solve the system of equations 

u'^ - xlli=o = ^Cfc,«(^,j)/^A;,Z (19) 

k,l 

where the coefficients Ck,i{i,j) are defined by 

^fij - rfij = ^Ck,iii,j)t^x\ 

k,l 

with A defined in ([5]) -([6]) with A replaced by A. 



3.2 Linear programs 

By optimizing over the pair of measures that satisfies the adjoint equations, the 
value V can be bounded, as follows: 

inf L{u, fi) < V < sup Li^v, fi) (20) 

which concerns linear programs over the measures, since L is the linear functional 
of the moments of ly and fj, given by 

ij m ij m 

and the infimum and the supremum are taken over the pairs of measures (z/, /i) 
supported on ([0, T] x R\i?, [0, T] x B) that satisfy the linear adjoint equations 
derived before. To formulate these optimization problems completely in terms of 
moment sequences, we need to express the condition that fj, and v be measures in 
terms of their moments, as in general there is no guarantee that any solution of the 
system (jl7p is the moment sequence of some measure. The problem to determine 
whether a given sequence is the moment sequence of some measure and, if so, 
whether this measure is uniquely determined (in which case the measure is called 
moment-determinate) has been extensively studied. It is known, see [19], that the 
Cramer condition 

e^^^^m{dx) < +00, for some c > 
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is a sufficient condition for a measure m to be moment determinate. In particular, 
any measure with compact support is moment-determinate. Further, the follow- 
ing Hausdorff conditions are necessary and sufficient for a given sequence to 
correspond to a moments of a measure m with support on the interval [a, b] (see 
e.g. 1): 

X^ r)(-iymj+fe >0 Vn,/c = 0,1,2,... (21) 
where the rhi are linear combinations of the mj, as follows: 



I 



mi = [0 — a 

i=0 



In fact, the m, are themselves the moments of a measure rh that is the affine 
transformation of m supported on [0, Ij. That these conditions are necessary 
immediately follows by observing that Jg y'^(l — y)""m{dy) is non-negative and by 
expressing rfij in terms of mj. More generally, given an array {mij,i, j = 0,1, . . .), 
the two dimensional Hausdorff-conditions (e.g. jl8j ) 

m n / \ / \ 

^ ^ / m \ / n \ > o_ Vn, m. A;, / = 0, 1, 2, ... (22) 

1=0 j=0 ^ * ^ ^"'^ 

where the are related to the niij by 

= {b- ar\d - cr' E C) i-^f-'i-cy-'m,,, 

i=Oj=0 V-^^ 

are necessary and sufficient conditions to guarantee that there exists a measure m 
supported on [a, b] x [c, d] such that rriij = J^^ x^y^m{dx, dy). See [18] for proofs 
and further background on problems of moments. 



3.2.1 Unbounded support 

In the case that the measure has unbounded support there also exists conditions 
to characterize a sequence of moments. These conditions are no longer linear but 
can be conveniently be formulated in terms of so-called moment and localizing 
matrices (their definitions are recalled in the appendix). For a sequence to be 
equal to the moments of some measure it is necessary and sufficient that these 
matrices are positive definite. See [12] or [6] for a proof of this fact. 



3.3 Approximations and convergence 

To be able to calculate lower and upper bounds for the value v we approximate the 
optimization problems in ()20|) by restricting the total number of moments used 
to N . If i3 is a finite interval, employing the moment conditions ()2ip and ()22p 
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results in the following (finite) linear programming problems: 



max 



mm 



4? 



i,j m 

subject to 



i.j m 



x^li=o = ^Cfc,/(«, j)/^fc,«, i + j<N, k + l<N 

(m) 



kJ 



with „ = Y^„("'\ /i 



conditions ([21])/ (1221) for u^j', /i. 



(m) (m) 



i + j<N 



In the case that the set i? is a half-line, the measures in question will not have 
bounded support and as a consequence in the above optimization problem the 
linear moment conditions (|2ip / (j22p are replaced by the quadratic moment con- 
ditions described in Section I3.2.11 The resulting optimization problems are then 
semi-definite programming problems. [13J provided convergence results for this 
SDP approach in a diffusion setting for Asian, European and single barrier op- 
tions. Restricting to the case that i? is a finite interval we show that the values 
v^\v^J^^ of the linear programs converge: 



Proposition 2 Suppose that the system (jl9p has a unique solution and that B is 
a finite interval. Then 



v_ \ V ana v , I v. 



as N 



oo. 



Remark. The presented approach can in principle be extended to a multi- 
dimensional jump-diffusion X with polynomial coefficients. For example, if B 
is a hyper-cube the adjoint equations and the moment conditions take analogous 
forms. The limitation in practice will be the capacity of the LP and SDP solvers 
to deal with large size programs. 



4 Numerical examples 

For the numerical examples we have used Matlab and the LP solver lp_solve. 
The problems were set up in Matlab and then solved using the Matlab interface 
to lp_solve. The numerical outcomes were compared with Monte Carlo results, 
implemented in Matlab using the Euler scheme. 

We will illustrate the method by valuing four different options. 

4.1 A double knockout barrier option driven by the Geometric 
Brownian motion. 

In this benchmark example we consider a European double knock-out call option 
with underlying St assumed to evolve as a geometric Brownian motion. The value 
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V of such an option is given by 



V = e"*^^ E[{St - i^)+/{^>T}] where 
T = mi{t>0:St^[Bd,Bu]}. 

For the ease of notation we will now drop the discounting e"''^. As a geometric 
Brownian motion has continuous paths, we know that the time-space process (t, St) 
will exit [0, T] x [B^i, By] either if S hits one of the barriers B^ or B^ or maturity 
is reached, so that the support for the exit location measure v is 

n = {[0,T] X {Bd}}U{[T] X [Bd,Bu]}U{[0,T] x {By}}. 

The set il. is partitioned into four parts with the restricted measures 

u^^^ and v^"^^ with support on [0, T] 

u^^^ and u^'^^ with support on [B^, K] and [K,Bu] respectively. 

The expected occupation measure /x is supported on the domain [0, T] x [Ba, B^] 
- See also Figure [U for an illustration. 







1^ 



(1) 



,(2) 



Bu 

(4) 



V 

K 



(3) 



V 

Bd 



Figure 1: Domain for the measures of double knockout option driven by the Geo- 
metric Brownian motion. 



Here the line-segments [B^, K] and [K, By] were chosen in such a way that the pay- 
off function restricted to each of those line segments is zero or linear. Further note 
that the measures v^^^ can all be characterized by Hausdorff moment conditions, 
since they are supported on line-segments. The value of the option is then given 

by 

v= f {x- K)+I{t=T}v{dt,dx) = vf' - Kv''^\ 
Jn 
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Case 1 


b = 0.1 


cr = 0.1 






Degree of moment 


9 


10 


11 


12 


Upper Bound 


0.9250 


0.9211 


0.9182 


0.9161 


Relative Error 


1.61% 


1.18% 


0.86% 


0.64% 


Lower Bound 


0.9096 


0.9100 


0.9102 


0.9103 


Relative Error 


U.UoTo 


U.UoTo 


n noo/ 
U.UzTo 


o m 0/ 


Exact solution 


0.9103 








Case 2 


b = 0.2 


cr = 0.2 






Degree of moment 


8 


9 


10 


11 


Upper Bound 


1.1656 


1.1611 


1.1569 


1.1534 


Relative Error 


2.06% 


1.66% 


1.29% 


0.99% 


Lower Bound 


1.1064 


1.1163 


1.1256 


1.1293 


Relative Error 


3.13% 


2.26% 


1.45% 


1.13% 


Exact solution 


1.1421 









Table 1: Numerical results for the Double knockout Barrier option with the un- 
derlying modelled by the Geometric Brownian motion. The option parameters are 
Bu = 5, Bd = l, K = 1.3, xo = 2, to = and T = 1. 

Using the form ([7]) of the infinitesimal generator of the Geometric Brownian mo- 
tion, the basic adjoint equation for this problem can be seen to be 

nfin-i,m - (^bm + Ym(m - 1)^ Hn,m = t^x^. 

This is valid for all n, m such that n + m < N , when we are using all moments 
up to degree N. To complete the setup of the problem we add the LP moment 
conditions for the measures u^^^ and fx with support as given above. 

The numerical results for two given sets of parameter values are given in Table 
[TJ We can see that we get fast convergence to the exact solution, which was 
calculated using the formula from |15j . 

4.2 A double knockout barrier option driven by the Variance 
gamma process. 

In this example we consider again a double knock-out option but now driven by 
a Variance Gamma process, which as described in Section [2j Since the Variance 
Gamma process is a finite activity jump process, it will not hit the barrier but 
jump across it. As a consequence the exit location measure v is supported on 

n = {[0,T] X [B„,oo)}U{{r} X [Bd,B^]}U{[0,T] X (-oo,5rf]}. 

In order to be able to calculate the value v of the option using the LP moment 
conditions, we will adjust the Levy measure rj as described in Section 13.11 to 
achieve bounded support. In this case we observe that any jump with absolute 
size larger than L = B^ — Bd will trigger an immediate knock-out. We note that 
the probability that no such a jump occurs before maturity is = e~^*"^ where 
A* = ry(M\[— -L, L]) with r] the Variance Gamma Levy measure given in ((HI). We 
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B +L 

u 



,0) 



B 



1^ 



V 

K 



u 

(4) 



V 

B 



(3) 
d 



,(2) 







T 



Figure 2: Domain for the truncated measures of a double knockout option in the 
Variance Gamma case. 

thus truncate the Levy measure rj by restricting it to absolute jump-sizes smaller 
L: 

f]{dx) = /||^|<i}r?(d3;). 
The value of the option in terms of moments of these measures is then 



V = p* / {x - Ky Ii^t^j,^v{dt, dx) = p^l^l 



(4) 



(4)i 



where the support of the truncated exit location measure is 

fi = {[0,r] X [Bu,Bu + L\}\j{[T] X [S^,SJ}U{[0,T] x [3^-1,3^]} 
and the four restrictions of v are 



(1) 

.(2) 



supported on [0, T] x Bu + L] 



v^'' supported on [0, T] x [Bd — L, B^] 
u^^^ and z^*-^-* supported on [B^, K] and [K, By] 

The domain of the truncated measures is shown in Figure [2l Denoting by 



c{k) = I ^y^k{y)dy 



(23) 



the moments of the truncated Levy measure fj and taking note of the form ([9]) the 
infinitesimal generator, we find the basic adjoint equation 



n,m— 1 



lyW + 1,(2) I / (3) , (4)\ _ ^„ ^ _ 

'^n,m ' '^n,m ' I '^m ' '^m J "■H'n-l,m UlllfJ,. 
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Case 1 


G = 8 


M = 12 


= 1.0000 




Degree of moment 


7 


8 


9 


10 


Upper Bound 


0.5045 


0.5030 


0.5022 


0.5017 


Lower Bound 


0.4946 


0.4983 


0.4987 


0.4994 


Relative Error 


0.13% 


0.09% 


0.05% 


0.07% 


Upu iime 


i.ZDZS 


O.Z4t)S 


b.ZODS 


i ( .aoos 


Monte Carlo 


0.5002 


Std Error 


0.0005 




Case 2 


G = 4 


M ^ 10 


= 1.0000 




Degree of moment 


6 


7 


8 


9 


Upper Bound 


0.5158 


0.5151 


0.5135 


0.5115 


Lower Bound 


0.4857 


0.4886 


0.4943 


0.4958 


Relative Error 


U.oa 70 


U.i I /O 


U.Z470 


n 1 no/ 


Monte Carlo 


0.5027 


Std Error 


0.0008 




Case 3 


G = 8 


M = 8 


= 1.0000 




Degree of moment 


5 


6 


7 


8 


Upper Bound 


0.5133 


0.5078 


0.5049 


0.5033 


Lower Bound 


0.4682 


0.4894 


0.4917 


0.4957 


Relative Error 


i. ( 1/0 


U.i47o 


n ono/ 
U.zUto 


U.U470 


Monte Carlo 


0.4993 


Std Error 


0.0006 




Case 4 


G = 3 


M = 6 


= 0.9998 




Degree of moment 


6 


7 


8 


9 


Upper Bound 


0.5277 


0.5237 


0.5197 


0.5182 


Lower Bound 


0.4672 


0.4720 


0.4745 


0.4772 


Relative Error 


1.32% 


1.24% 


1.39% 


1.27% 


Monte Carlo 


0.5041 


Std Error 


0.0011 





Table 2: Numerical results for the Double knockout barrier option driven by the 
Variance Gamma process. The option parameters are in all cases, Bu = 1, B^i = 
— 1 and K = —0.3, the model parameters are b = 0.2, xq = 0, T = 1 and C = 0.5. 
The truncation size is L = 2 in all cases, and p^, is the probability that no jump 
of absolute size larger than 2 occurs before T. 

for all n, m such that m + n < N . As before, to complete the LP problem we add 
the appropriate LP moment conditions for the measures i/^*^ and fi with support 
as given above. 

Numerical results can be found in the Table [2l The relative error was calculated 
using the arithmetic mean of the upper and lower bounds and the Monte Carlo 
outcome (taking the latter as the 'true' result). Studying the results we see that we 
get tight bounds within 8 or 9 moments. Beyond 10 or 11 moments we experienced 
instabilities with the LP solver. 

In Case 1 the execution times are shown, which should be compared to the ex- 
ecution time for the Monte Carlo simulation that was around 52 minutes. As 
we employed a basic Euler scheme for the Monte Carlo simulation the speed of 
convergence of the Monte Carlo simulation could be improved by using more spe- 
cialised Monte Carlo schemes and also by changing to a compiling programming 
language. However considering the time difference, it should be clear that the 
method of moments will still be considerably faster. 
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4.3 The American Corridor under a CIR model with constant 
interest rate 



An American corridor is a contract traded in the Foreign exchange markets that 
pays a continuous rate until either the underlying leaves the corridor or maturity 
is reached, whichever comes earlier (see for example |20j or [21] for background). 
The value of this contract is given by 



E 



/ e ^^dt = / fi{dt, ds) = fiofi where 
Jo \ Je 



T = T Amf{t>0: St^ [Bd,Bu]}, 

with 6 = {[0, T] X [BdjBy]}- We will model the underlying as a Cox Ingersoll 
Ross (CIR) process, evolving according to the SDE (fTO]l . Since the CIR model is 
continuous the supports of the different measures are given as follows: 

u^^^ and u^'^^ supported on [0, T] 

u^^^ supported on [B^, By] 

H is supported on [0,T] x [Bii,Bu] 

In view of the form of the infinitesimal generator for the CIR process pi|) we can 
now assemble the basic adjoint equation for this problem, 



m 

2 



{am + r)fXn,m - ( abm + ^m{m - 1) ) fj,n,m-i = t^x^ 



for all m, n such that m + n < N , and add appropriate LP moment conditions as 
before. 

The results are reported in Table [3l We observe that tight bounds are achieved, 
and that in cases 1 and 2 the upper bounds are accurate for 9-10 moments, with 
relative errors 0.16% and 0.12%. We also see that the speed of convergence varies 
with the particular parameter values. 



4.4 Double No Touch option under the exponential Variance Gamma 
process with a non constant interest rate 

A double no touch option pays one unit at maturity if the underlying has not 
crossed either of the barriers B^ or B^. Its value can be expressed as 

V = i?[e~"^/{T->r}] where 
T = mi{t>0:St^[Bd,Bu]}. 

Letting the underlying St be an exponential Variance Gamma process we note 
that the stopping time r is equivalent to 

r = inf{t >0:Xt^ [log(i?rf), log(S„)]} 
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Case 1 


(7 = 0.2 


r = 0.1 






Degree of moment 


10 


11 


12 


13 


Upper Bound 


0.9516 


0.9516 


0.9516 


0.9516 


Lower Bound 


0.9274 


0.9345 


0.9391 


0.9421 


Relative Error 


1 1 00/ 


n V/i 0/ 
U. / 4/0 


U.oU/o 


Q /I 0/ 


Monte Carlo 


0.9501 


Std Error 


0.0002 




Case 2 


a = 0.2 


r = 0.05 






Degree of moment 


9 


10 


11 


12 


Upper Bound 


0.9754 


0.9754 


0.9754 


0.9754 


Lower Bound 


0.9394 


0.9504 


0.9577 


0.9624 


Relative Error 


1 TOO/ 

1.737c 


1.1d7o 


U.Ty/o 


n c /1 0/ 
0.5470 


Monte Carlo 


0.9742 


Std Error 


0.0002 




Case 3 


a = 0.3 


r = 0.1 






Degree of moment 


11 


12 


13 


14 


Upper Bound 


0.9343 


0.9325 


0.9315 


0.9307 


Lower Bound 


0.8961 


0.9024 


0.9067 


0.9095 


Relative Error 


0.76% 


0.52% 


0.34% 


0.29% 


Monte Carlo 


0.9222 


Std Error 


0.0011 





Table 3: Numerical results for the American Corridor modelled by the Cox Inger- 
soll Ross model. The problem parameters are = 0.5, = 1.5, xq = 1 and 
T = 1, for the model parameters a = 0.5 and 6=1 are fixed. 

where Xt is a Variance Gamma process. Under no-arbitrage pricing the process 
e'^^^St = e~°'*~^^* needs to be a martingale which is equivalent to the requirement 
that 

/oo 
(e^ - l)r]{dx) = r{t) 
-oo 

so that b(t) is determined by our choice of r{t). For the Variance Gamma process 

c = ly - m^) = c (log (^) + log (^)) 

With these points in mind we find, 

{A - r{t))f{t, x) = ^ + b{t) + j[f{t,x + y)- f{t, xMdx) - r{t)f{t, x) 

= % + {r{t) - + / ^ + y) - ^Mdx) - r{t)f{t, x) 
We have chosen to study interest rates of the type 

r{t) = n + Vst^. 

Since there is no need too split the exit location measure at maturity, the domain of 
u only needs to be split into three parts Ki = [0, T] x [0, B^], K2 = [0, T] x 00) 
and = {T} x [Bd^Bu]., so that in this case the value v of the option can be 
expressed in terms of moments as follows 

V = E[e-'^^I{T >T)]= [ I{t = T)v{dt, ds) 

Jn 

JKa 
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Case 1 


n = 0.05 


rs = 0.05 






Degree of moment 


6 


7 


8 


9 


Upper Bound 


0.9356 


0.9355 


0.9355 


0.9355 


Lower Bomid 


0.8453 


0.8757 


0.9042 


0.9143 


Relative Error 


4. I\)/q 


Q 1 70/ 
d.i I /O 


1.04/0 


1 1 no/ 
i.iU/o 


Monte Carlo 


0.9352 


Std Error 


0.0002 




Case 2 


n = 0.05 


rs = 0.1 






Degree of moment 


6 


7 


8 


9 


Upper Bound 


0.9203 


0.9201 


0.9200 


0.9200 


Lower Bound 


0.8196 


0.8533 


0.8836 


0.8957 


Relative Error 


5.O<370 


3.0070 


1 m oy 
1.9l7o 


i.Z67o 


Monte Carlo 


0.9194 


Std Error 


0.0002 




Case 3 


n = 0.1 


rs = 0.1 






Degree of moment 


7 


8 


9 


10 


Upper Bound 


0.8752 


0.8752 


0.8752 


0.8752 


Lower Bound 


0.7980 


0.8319 


0.8449 


0.8565 


Relative Error 


6.68% 


4.73% 


2.84% 


2.09% 


Monte Carlo 


0.8746 


Std Error 


0.0002 





Table 4: Numerical results for the double no touch option with barriers at = 2 
and Bii = 0.5 and maturity T = 1, driven by an exponential Variance Gamma 
process with parameters C = 0.5, G = 8 and M = 12, and with Sq = 1. Note that 
the upper bounds are very accurate already for 6-7 moments, with relative errors 
less than 0.1%. 

where as before i/*^*-* = u{- D Ki). In terms of moments the basic adjoint equation 
is then given by 

-1-7/(2) A_T'^,X3) _ m _ —ri.li — r ii ,o 

'^^n,m ' '^n,m ' ^ '^m ''O'^O — "-f-n— l,m ' bh'"n,m ' sh''n+2,'m 

+ {rb - c)mfXn,m-l + rsmfXn+2,m~l + X] ( t ]c(fc)^n,m-fc, 

k=l ^ ^ 

where c{k) is given in ()23p . 

The numerical results are presented in Table HI In all cases we see tight bounds 
nicely agreeing with the Monte Carlo simulation result. We also observe that the 
upper bound is very accurate already for a small number of moments. For example, 
for 7 moments, the relative errors of the upper bounds in the three different cases 
are 0.032%, 0.076% and 0.069%, respectively. 



5 Conclusion 

We have presented a method of moments approach that can be used to price double 
barrier-type options driven by 'polynomial' jump-diffusions, allowing for a non- 
constant (deterministic or stochastic) interest rate. An infinite-dimensional linear 
program was derived, which was then approximated, depending on the choice 
of moment conditions, either by a sequence of LP problems, or by a sequence 
of SDP problems. Although the SDP-type problems may be theoretically more 
appealing as the SDP method naturally handles measures with unbounded sup- 
port, further development of stable SDP solvers would be needed for this method 
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to be truly usable in practice. Since, on the other hand, the LP solvers are in 
a more advanced state of development and several (commercial) LP solvers are 
available capable of solving (large scale) LP problems, we focussed on the LP 
approach. We formulated the approximating programs as LP problems by using 
truncation, and provided theoretical convergence results for this approach. We 
illustrated the method with numerical examples, using the Matlab interface of the 
solver lp_solve, and compared the outcomes with Monte Carlo simulation results. 
We found that accurate results with tight upper and lower bounds were obtained 
with a small number of moments in most of the examples, and observed that the 
algorithm was significantly faster than Monte Carlo simulation. 
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A Proofs 



A.l Proof of Proposition [2] 



Since the number of equations grows with A^, it follows that v^^^ is monotone 
increasing, since the minimum taken over a smaller set of elements is larger. 

Note that v is a finite linear combination of moments. Because of the fact that 
the support of the different measures is compact, if follows that each moment is 
bounded, so that v is bounded and is the minimization of a linear function 

over a bounded set. Thus, the minimum is finite and attained at a vector 
= {qj^)i that satisfies the corresponding linear system of equations. Thus, for 
each fixed i there exists a such that, for along a subsequence, qf — > q^. In 
fact, by a diagonal argument it follows that there exists a subsequence such 
that, as ^ oo, 



q* for all i. 



Clearly, q* satisfies the infinite system and thus under the assumption that there 
exists a unique sequence that solves the infinite system, it follows that q* must 
be equal to Moreover, ^ and v are the unique measures corresponding 

to these moments, as they are both moment-determinate. Thus, v = L{q*) and 

V^!"^ T L{q*) 

The proof of the convergence of the sequence {v^^'^ ] 



is similar and omitted. 



□ 



A. 2 Proof of Proposition [T] 

We will show the following lemma: 

Lemma 1 For any bounded stopping time r and f G C^'"^ with 



E 
with 



a 



•.dl 
dx 



it,Xt 



+ E 



e~''^\g{t,Xt,y)\A{dy)dt 



< +00, (24) 



df 

9{t, X, y) = f{t, X + yX{t, x)) - f{t, x) - —{t, x)X{t, x)yl\y\^i, 
eqn. (fT5|l holds true. 



The proposition is a direct consequence of this lemma, since under the condition 
(|18p the integrability conditions are satisfied for each monomial t^x^ . 

Proof of Lemma Ul' Applying (a general form of) Ito's lemma to the stochastic 
process e~°'^ f{t, Xt) (which is justified as / S C^'^) shows that 

e-''^f{t,Xt)-f{0,Xo) = Mt+ f e'^'iAf -rf){s,Xs)ds (25) 

Jo 

where Af is given in @ and Mt is the local martingale given by 

Mt = f e-"V(t, Xt)dWt + f e-"^g{t, Xt,y)(P{dy, dt), 

Jo J[0,t]xR 
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where cj) denotes the compensated jump measure associated to J (with compen- 
sator K{dy)dt). Note that the identity (j25p remains vahd with t replaced by t A r 
with M"^ = {Mtfsr^t > 0} a local martingale. Under the integrability conditions 
(I24|) it holds that that M"^ is a zero mean martingale. Taking expectations in (1250 
shows thus that 



E[e-^^f{t,Xt)]-f{Q,x) 
which shows that (jlSp is valid. 



E 



e-'^^iAf -rf){s,Xs)ds 



□ 



A. 3 Proof of Corollary [T] 



Since the jump-diffusion X with drift 6, volatility o", Levy measure A and discount- 
ing f satisfies (llSp , Proposition [T] yields that the measures v and corresponding 
to X satisfy 



oli=o = '^Ck,i{^,j)J^Ki i,j = 0, 1,2, 



kd 



where 



{Afij - rfij){t,x) =^Ck,iii,j)t^x'', 



k,l 



with fij = fx^. To complete the proof we will now show that 

f = y h{t,x)i>{dt,dx) + J g{t,x)fi{dt,dx). (26) 

Denoting by p the first time that a jump of J of size smaller than L_ or larger 
than L_|_ and let 

f = inf{t >0:Xt^B}. 

Then, if p > r it holds that Xt/\T = Xt/sr for all t > and in particular 

T = f and Xj- = Xj- = Xf. 

Also, since it is assumed that h{t,x) = for x ^ B, we have that /i(r, X,-) = 
h{T, Xr)l^pyr}- Taking note of these observations, it follows that 

E[e-''^h{T,Xr)] = E[e-"-h{T,Xr)l{p>r}] 
= E[e-"^h{f,Xf)^p^,}] 

= E[e-''^''^^h{f,Xf)] 
= E[e-''^h{f,Xf)] 

where A := A(M\[L_, and we used that p follows an exponential distribution 
with mean A~^, independent of X. Similarly, 



E 



I e"^g{s,Xs-)ds = E [ e-"=g(s, X,_)l|,<^|ds 
Jo J Uo 



E 
E 
E 
E 



e °''g{s,Xs-)l{s<p}ds 
e-'''-~^'g{s,Xs-)ds 
e-''^gis,Xs-)ds 
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The two identities imply that ()26p holds true, and the proof is complete. 



□ 



B Semi-definite moment conditions 



The moment matrices are defined as follows: 



Moment Matrices Let 



(x", |q| < k) — (1, Xi, . . . , Xn, x\, X1X2, ■ ■ ■ , Xi, ^X2, ■ ■ ■ , X^) 



(27) 



be the usual basis of polynomials in n variables with degree at most k. 

Given a series of moments of a measure m = {ma, a £ N"} let m = {mi,i G N} 
be that sequence ordered in accordance with (|27p . The Moments matrix Mk{m) 
is then defined as 

Mfe(m)(l,i) = Mk{m){i, 1) = mj_i,for i = 1,. . . ,k + 1, 
Mfc(m)(l, j) = nia and Mk{m){i, 1) = m/3 Mk{m) = rua+p 

where Mk{m)(i, j) is the (z, j)-entry of the matrix Mk{m). 

Localising Matrices Given a polynomial q with coefficients (qa) in the basis 
()27p . If (3{i,j) is the /5 subscript of the (z,j)-entry of the moment matrix Mk{m) 
then the localising matrix is defined by 



The choice of the function q depends on the support of the measure. For example, 
in the one-dimensional setting we have three cases, 

1. Support on [a, 6], with function q = {h — x){x — a) 

2. Support on [a, cxd), with function q = {x — a) 

3. Support on {—00, a], with function q= [a — x) 

In terms of moment and localizing matrices the characterization is then as follows: 
Given m = (mo, mi, . . . , m2r) the condition that Mr{m) and Mr-i{q, m) are pos- 
itive semi-definite are sufficient conditions for the elements of m to be the first 
2r + 1 moments of a measure supported on the appropriate interval. 




a 
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